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Summary 


Ecological niche modeling has allowed several advances in the process of species delimitation. In the present study, 
I used this method to evaluate the climatic divergence between two scorpion species, Mesobuthus eupeus and M. 
phillipsii. The ecological niche models (ENMs) were created based on presence-only data with the maximum 
entropy method. The created models, results of spatial analysis (PCA and Hotelling discriminant), and an identity 
test suggested that the divergence between these two species is associated with significant divergence in their 
ecological niches. The results of this study provide additional support for the taxonomic validity the studied species. 





Introduction 


The species concept and the process of species 
delimitation have been a matter of controversy since the 
early days of systematic biology, and several new 
methods have been proposed, developed and compared 
to accelerate the process (Wiens & Penkrot, 2002). 
Among these methods, ecological factors should help 
delimit species, assuming each species has developed 
and occupied its own ecological niche (Kozak et al., 
2008; Rissler & Apodaca, 2007; Kozak & Wiens, 2006; 
Phillips et al., 2006; Wiens, 2004). 

The suite of methods variously called species 
distribution modeling, habitat modeling or ecological 
niche modeling (ENM) have been frequently applied to 
identify places suitable for the survival of populations of 
a species through identification of their environmental 
requirements (Soberón & Nakamura, 2009). ENM is a 
rich area of study that has seen tremendous growth in 
past years and is essential to diverse applications in 
ecology (DeVaney, 2008; Mingyang et al., 2008), 
conservation biology (Papes & Gaubert, 2007; Urbina- 
Cardona & Flores-Villa, 2009; Thorn et al., 2009), 
vector and disease control (Ayala et al., 2009; Chamaillé 
et al., 2010; Foley et al., 2010; Saupe et al., 2011), and 
evolution and systematics (Rissler & Apodaca, 2007; 
Raxworthy et al., 2007; Makowsky et al., 2010). These 
methods have allowed several advances in under- 
standing the speciation process and geographic ecology 
of species, as well as ecological and evolutionary 


determinants of spatial patterns of biological diversity 
(Peterson & Vieglais, 2001; Stockwell et al., 2006; Elith 
et al., 2006). Although the systematic application of 
ENM is well known, there has been little discussion of 
its application on species delimitation (Elith et al., 
2006). A variety of methods exist for modeling the 
potential distribution of species based on occurrence 
data and digital maps of environmental variables (Elith 
et al., 2006). These models employ relationships be- 
tween environmental variables and known species' 
occurrence localities to describe abiotic conditions 
within which populations can be maintained (Kozak & 
Wiens, 2006; Raxworthy et al., 2007). 

The species Mesobuthus eupeus (C. L. Koch, 1839) 
occurs in the Palearctic region from eastern Anatolia to 
China (Fet, 1994; Karataş & Karataş, 2001, 2003; 
Teruel, 2000; Gromov, 2001; Qi et al., 2004; Shi et al., 
2007). This species represents at least 14 morphological 
subspecies that can be distinguished on the basis of 
coloration, color patterns and carination (Fet & Lowe, 
2000; Mirshamsi et al., 2010). Nevertheless, con- 
siderable debate exists regarding the nature of these 
morphological variations (Fet, 1994): Do the various M. 
eupeus subspecies simply represent different morphs of 
a single widely distributed species? Do the subspecies 
represent distinct evolutionary entities that could be 
considered as a separate species? 

The results of recent studies based on mitochondrial 
sequence data and multivariate statistical analyses 
(Mirshamsi, 2010; Mirshamsi et al., 2010) clearly re- 


vealed two divergent lineages within M. eupeus (~6% 
divergence between COI sequences). The geographic 
distributions of these two lineages are consistent with 
the proposal that the formation of Zagros Mountains in 
the south, and consequent uplifting of Alborz Mountains 
in the north, may have influenced their evolution. There- 
fore, this provided support for raising M. e. phillipsii 
(Pocock, 1889) to species level (Mirshamsi et al., 2010; 
Mirshamsi, 2010). The results of a detailed morpho- 
logical analysis of this species also strongly support 
elevating the southwestern subspecies to species status 
(Mirshamsi et al., 2011). 

In this study, I examined the association of climatic 
variation and geographic isolation between M. eupeus 
and M. phillipsii. I combined presence-only occurrence 
data, digital maps of climatic variation, and a maximum 
entropy method to estimate the potential distributions of 
these species in regard to their taxonomic status. Using 
this method, I also tested the hypothesis that historical 
and current geographic isolation between these species is 
a result of niche divergence. 


Material and Methods 


Occurrence data 


Distributional data for M. eupeus and M. phillipsii 
were obtained in the field and from literature records 
that could be georeferenced (Pocock, 1899; Birula, 
1903, 1905; Navidpour et al., 2008 a, b, c, d; Pirali- 
Khairabadi et al., 2009; Navidpour et al., 2010; Navid- 
pour et al., 2011; Kovařík et al. 2011). Literature 
records were evaluated to check the species ident- 
ification and localities, and all localities with uncertainty 
greater than 5 km were excluded from the analysis. 


Ecological niche modeling and data analysis 


Niche modeling requires two types of data: coord- 
inates (georeferenced occurrence records) and GIS- 
based maps of the environmental variables (e.g. climatic 
variables such as temperature and precipitation) that are 
likely to affect the suitability of the environment for a 
given species (Kozak et al., 2008). I used the program 
Maxent ver. 3.3.2 to model the potential distribution of 
M. eupeus and M. phillipsii. Maxent combines occur- 
rence records, random background points and environ- 
mental layers to provide a gridded model of potential 
distribution of the studied species using a statistical 
approach known as Maximum Entropy. Maxent was run 
using a subset of bioclimatic layers (out of 19) repre- 
senting trends in temperature, precipitation, and 
seasonality with 30 arc second resolution (~1 km”) 
extracted from the WorldClim database (Hijmans et al., 
2005; http://www.worldclim.org). 
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To avoid over-fitting problems which may result 
from inclusion of too many climatic variables, the 
environmental information from 1000 randomly gen- 
erated geographic points across the study area was 
extracted using ArcMap. Correlation matrices were then 
created for all 19 climatic variables within each of the 
temperature and precipitation categories. A Pearson 
correlation coefficient of 0.75 (using PAST software, 
Hammer et al., 2001) showed that the variables were 
highly correlated (Rissler et al., 2006). For such 
correlated variables the more biologically informative 
variable was used to create niche models (Rissler & 
Apodaca, 2007). Finally, six temperature and five 
precipitation variables were selected and used to 
generate ENMs. The selected variables included: BIO1= 
annual mean temperature; BIO2 = mean diurnal temp- 
erature range; BIO3 = isothermality; BIO7 = annual 
range in temperature; BIO8 = mean temperature of 
wettest quarter of the year; BIO9 = mean temperature of 
driest quarter of the year; BIO15= precipitation of 
seasonality; BIO16= precipitation of wettest quarter of 
the year; BIO17= precipitation of driest quarter of the 
year; BIO18 = precipitation of warmest quarter of the 
year; BIO19= precipitation of coldest quarter of the 
year. 

Models for targeted species were created using 
a total of 365 point localities. All Maxent runs were 
adjusted with a convergence threshold of 1.0E-5 
with 1000 iterations and the regularization value 
was set to 0.1. I also chose the logistic output 
format, displaying suitability values from 0 (un- 
suitable) to 1 (optimal). Finally, I imported the 
resulting ASCII files into DIVA-GIS to visualize 
the models. 

Evaluation of predictive models was conducted 
by splitting localities into 75% training and 25% 
test data partitions, and five replicates of random 
data partitions were performed (Phillips et al., 
2006). Model validation was performed by calculating 
the area under the curve value (AUC), which reflects the 
model's ability to discriminate between presence records 
and random background points. AUC values range from 
0.5 for models without any predictive ability to 1.0 for 
models with perfect predictive ability. Therefore, AUC 
values >0.9 considered to have 'very good', >0.8 'good' 
and >0.7 'useful' discrimination ability (Hawlitschek et 
al., 2011). 

Additional statistical analysis of modeling results 
was performed using ENMtools (Warren et al., 2010). 
Niche overlap of M. eupeus from Central Asia and 
eastern and northern Iran and of M. phillipsii from Iraq, 
Syria, southeast Turkey and southwest Iran was 
calculated using Schoener's D and the J statistics. Both 
of these indices range between 0 and 1, and values close 
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Figure 1: Ecological niche models constructed using Maxent. (A) M. eupeus niche model (AUC=0.966+0.012); the resulted 
ENM includes the type localities of five morphological subspecies, M. e. eupeus, M. e. philippovitschi, M. e. afghanus, M. e. 
thersites and M. e. haarlovi. (B) M. phillipsii niche model (AUC=0.977+0.013); the resulted ENM also includes the type 
localities of morphological subspecies M. e. kirmanensis and M. e. mesopotamicus. White circles indicate the occurrence points 
of each species. 
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Figure 2: Result of the identity test performed by ENMtools. The red and black arrows indicate the result of ENMtools’ niche 
overlap analysis representing the true calculated niche overlap (D and I). The red and black curves represent the niche overlap 
values created in the replicates of the identity test. The calculated overlap values (D and I) are outside the 99.9% confidence 
intervals of the identity test results and thus statistically significant. 


to 1 explain high similarity between ENMs (Schoener, 
1968). Also, the identity test included in ENMtools was 
performed. The niche identity test was used to determine 
whether the niche models generated for two species are 
identical or exhibit a statistically significant difference. 

Finally, to evaluate the overall level of divergence 
in the ecological niches of M. eupeus and M. phillipsii, a 
principal components analysis (PCA) was conducted 
using the extracted values for each species. I then used a 
multivariate Discriminant/Hotelling test with PCA axis 
scores as dependent variables and lineage as the fixed 
factors to determine whether the observed separation in 
the ecological niche was statistically significant. 


Results 


In total, 176 and 189 point localities were used for 
generating ENMs of M. eupeus and M. phillipsii, 
respectively. Ecological niche models for M. eupeus and 
M. phillipsii are presented in Figs 1A & B. According to 
the estimated average AUC values, 0.966 + 0.012 for M. 
eupeus and 0.977 +0.013 for M. phillipsii, the ability to 
differentiate presence from random background points 
for all generated models was larger than 0.9 and thus 
considered as "very good". 








A heuristic estimate of relative contribution of the 
climatic factors to the MaxEnt models showed that for 
the distribution of M. eupeus, “mean temperature of 
wettest quarter of the year”, “precipitation of coldest 
quarter of the year” and “precipitation of warmest quar- 
ter of the year” were the variables of highest importance 
with 27.9%, 25.5% and 12% of contribution, respect- 
ively. For M. phillipsii, “precipitation of warmest quarter 
of the year” was the variable of highest importance with 
34.8% of contribution. For this species, “precipitation of 
coldest quarter of the year” and “mean temperature of 
wettest quarter of the year” was the second and third 
most important predictors with 31.2% and 9.3% of 
contribution, respectively. 

According to the calculated niche overlap between 
M. eupeus and M. phillipsii, I=0.5783 and D=0.3185, the 
overlap between niche models is considered low. 
According to the result of identity test, which combines 
the occurrence data from two species into a common 
pool, the null hypothesis of niche identity was rejected. 
The results suggested that the estimated niche models 
for M. eupeus and M. phillipsii based on the climatic 
variables are significantly distinct (Dyj=0.8373+0.002 
SE vs. Dy=0.3185 and J p~o=0.9641+0.0044 vs. I 
H/=0.5783) (Fig. 2). 
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The result of PCA based on 11 climatic variables 
indicated that the two species were distributed in diver- 
gent environmental spaces (Fig. 3). The first three com- 
ponents resulting from this analysis explained 40.95%, 
21.87% and 17.82% of the total variation. The climatic 
variables BIO1, BIO2, BIO3, BIO8, BIO9, BIO15 and 
BIO19 have positive loadings on the first component 
while the climatic variables BIO2, BIO3 and BIO7 have 
positive loadings on the second component. The result of 
the Discriminant/Hotelling analysis based on PCA axis 
scores also revealed that the separation between ENMs 
was statistically significant (P«0.005). 

According to the model, climate suitable for M. 
eupeus ranges from Central Anatolia in the west to the 
east of Kazakhstan (Fig. 1A). The geographic distri- 
bution of this species appears to be limited between 
28°N and 46°N. The distributional model includes the 
terra typica of the nominotypical subspecies M. e. 
eupeus, as well as M. e. afghanus (Pocock, 1889), M. e. 
haarlovi Vachon, 1958, M. e. thersites (C. L. Koch, 
1839), and M. e. philippovitschi (Birula, 1905), 

Moreover, on the basis of the model for M. 
phillipsii, suitable climate conditions for this species are 
distributed from southeast of Turkey to southeast of Iran 
(Fig. 1B). The distributional model includes the terra 
typica of M. phillipsii, as well as two other morpho- 
logical subspecies, traditionally classified under M. 
eupeus: M. e. kirmanensis (Birula, 1900), and M. e. 
mesopotamicus (Penther, 1912). 


Discussion 


Species are essential units of study in biosystematic 
and evolutionary biology, and the precise delimitation of 
species, the process by which the boundaries of species 
are determined, is an increasingly important task of 
biosystematics (Wiens & Servedio, 2000; Wiens, 2007). 
In contrast to the major progress in the phylogenetic 
reconstruction process, there has been relatively little 
progress in the methodology of species delimitation 
process (Wiens, 2007). The use of ecological data and 
GIS-based analyses of climatic variables in determining 
the species boundaries, proposed by Raxworthy et al. 
(2007) and Rissler & Apodaca (2007), was a significant 
advance in species delimitation. 

Traditionally, M. eupeus was known as one of the 
most widespread buthid scorpions in the Palearctic 
region. The wide geographic distribution of this species 
is accompanied by a large amount of morphological 
variation, which has been used to delimit a dozen 
subspecies (Farzanpay, 1987; Fet, 1994). Since these 
subspecies have mainly been described on the basis of 
inconclusive external characteristics, the traditional 
taxonomy of this taxon had been complicated and some 
authors believed that M. eupeus is a species complex 
(Fet et al., 2003). 


Recently, the result of a molecular phylogenetic 
study of the M. eupeus species complex based on COI 
sequence data clearly indicated two distinctly divergent 
Northern (M. eupeus) and Southern (M. phillipsii) 
clades. The genetic divergence between these two line- 
ages was nearly comparable to that observed between 
congeneric species [i.e. 6% mtDNA sequence diver- 
gence (Mirshamsi et al., 2010)]. These two clades also 
showed different morphometric features based on uni- 
and multivariate statistical analyses. Therefore, in order 
to decide whether the two clades warrant separate 
taxonomic designations, a detailed revision of diagnostic 
characteristics of this taxon and its morphological sub- 
species was recommended (Mirshamsi et al., 2010). 
Finally, the morphological and morphometric exami- 
nation of different subspecies of the nominal M. eupeus 
clearly revealed that there are at least two valid species 
belonging to the M. eupeus species complex, including 
M. eupeus and M. phillipsii (Mirshamsi et al., 2011). 

ENMs based on 365 occurrence records and 11 
climatic variables clearly indicate that these two species 
differ in their ecological niches. Also, according to the 
results of the identity test, niche divergence between the 
studied species proved to be significant (Fig. 2). There- 
fore, the climatic envelope of M. eupeus differs from 
that of M. phillipsii. Moreover, the analysis of the spatial 
data (PCA) find that, in the genus Mesobuthus, closely 
related species occupy divergent ecological niches that 
show relatively little spatial overlap when projected in 
geographic space (Fig. 3). 

Our results suggest that M. eupeus and M. phillipsii 
are parapatric species with a contact zone formed at the 
border of their distributional ranges, which corresponds 
to the Zagros Mountain system. According to the ENMs, 
suitable climate for M. eupeus mainly occurs in arid and 
semi-arid regions of Iranian Plateau, Central Asia, and 
Central Anatolia. Notably, there is no confirmed record 
of M. eupeus from Central Iranian and Lut deserts, 
which is also predicted as unsuitable in the MaxEnt 
model of M. eupeus. The model clearly shows that the 
western border of the central Iranian desert and eastern 
foothills of the Zagros Mountains may act as geo- 
graphical barriers to dispersal of this species. According 
to the ENM and the type localities of traditional valid 
morphological subspecies, M. eupeus includes five 
morphological subspecies; M. e. eupeus, M. e. afghanus, 
M. e. haarlovi, M. e. philippovitschi and M. e. thersites. 
The subspecies M. e. philippovitschi has been synony- 
mized with the nominotypical subspecies, M. e. eupeus, 
by Kovarik et al. (2011) (Fig. 1). 

The data presented here also indicate that M. 
phillipsii mainly occurs in western foothills of Zagros 
Mountains towards the Tigris-Euphrate River drainage 
in Eastern Iraq. The model for M. phillipsii suggests that 
the niche of this species extends from Baluchistan in the 
southeastern Iranian Plateau, and westward to the Tigris- 
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Euphrates River drainage in eastern Iraq. For this 
species, the southern and western foothills of Zagros 
Mountains might act as the main geographic barrier to 
its dispersal. Notably, according to the ENM, the range 
of M. phillipsii includes also the type localities of 
subspecies M. e. kirmanensis and M. e. mesopotamicus 
which were traditionally classified under the name of M. 
eupeus. The results of present study suggest that these 
two subspecies should be transferred to M. phillipsii 
(Fig. 1). The subspecies M. e. mesopotamicus has been 
synonymized with M. phillipsii (under M. eupeus 
phillipsii) by Kovarik et al. (2011) who did not find 
considerable difference between topotypical populations. 

The ENMs predicted by the maximum entropy 
method and presence-only occurrence data represents 
the realized niche, because correlations among the 
species distribution and the environmental factors are 
obtained through the presence data (Sillero, 2011). The 
actual geographic distribution is a modification of fun- 
damental niche and could be defined according to the 
complex interaction of the realized environment, his- 
torical and biotic factors, as well as the fundamental 
niche. Because in presence-only niche modeling, we 
look for all the habitats suitable for the species, we 
cannot distinguish between really occupied and unoccu- 
pied habitats as the true absence data are not included in 
the models. This is perhaps the reason for the existence 
of false-positive predicted areas in the resulting geo- 
graphic distribution of both species. The results of 
analyses of spatial data all support the conclusion that 
the genetic and morphological divergences between the 
studied species are associated with significant diver- 
gence in the ecological niche of M. eupeus and M. 
phillipsii. 

According to the results of the current study, the 
present-day distribution of M. eupeus and M. phillipsii 
may have been affected by topographic barriers such as 
the Zagros orogeny. A model of vicariant speciation by 
the Zagros orogeny is consistent with the divergence 
time estimated based on mtDNA sequences (i.e. 4.60- 
3.15 mya) (Mirshamsi et al., 2010) and the distribution 
models of these species presented here. Based on this 
model, the evolution and speciation of these species and 
other terrestrial faunas in the Iranian Plateau have been 
influenced by uplifting of Zagros Mountains, which 
occurred in the late Tertiary approximately 10-5 mya 
(Macey et al., 1998; Gòk et al., 2003). In mid-Pliocene 
(5-3 mya), the Zagros orogeny was intensified as the 
Arabian plate rifted from Africa, leading to sinking and 
formation of the Central Iranian and Lut deserts. These 
geologic processes may have acted as the main vicariant 
events that partitioned Mesobuthus into two species; M. 
eupeus inhabiting northern regions of the Iranian 
Plateau, and M. phillipsii which became adapted to the 
southern regions of the Iranian Plateau, western foothills 


of the Zagros mountains and the Tigris-Euphrates River 
drainage in eastern Iraq. 
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